#include "borwein.h"

/* Armazena os valores iniciais para o método de Borwein */
void start_values_borw(mpf_t a, mpf_t y)
{
    /* a0 = 6-(4*sqrt(2)) */
    mpf_init_set_ui(a, 2); 
	mpf_sqrt(a, a); 
	mpf_mul_ui(a, a, 4); 
	mpf_ui_sub(a, 6, a);
    
    /* y0 = sqrt(2)-1 */
    mpf_init_set_ui(y, 2); 
	mpf_sqrt(y, y); 
	mpf_sub_ui(y, y, 1);
}

/* Cálculo de y */
void calc_y(mpf_t y1, mpf_t y0)
{
    mpf_t aux1, aux2;
    mpf_inits(aux1, aux2, NULL);

    /* aux1 = (1-y0^4)^(1/4) */
    mpf_pow_ui(aux1, y0, 4); 
	mpf_ui_sub(aux1, 1, aux1); 
	mpf_sqrt(aux1, aux1); 
	mpf_sqrt(aux1, aux1);
	
    /* aux2 = aux1+1 */
    mpf_add_ui(aux2, aux1, 1);
	
    /* y1 = (1-aux1)/aux2 */
    mpf_ui_sub(y1, 1, aux1); 
	mpf_div(y1, y1, aux2);

    mpf_clears(aux1, aux2, NULL);
}

/* Cálculo de a */
void calc_a(mpf_t a1, mpf_t a0, mpf_t y1, int k)
{
    mpf_t aux1, aux2;
    mpf_inits(aux1, aux2, NULL);

    /* Calcula aux1 = y1*(1+y1+y1^2) */
    mpf_pow_ui(aux1, y1, 2);
	mpf_add(aux1, aux1, y1); 
	mpf_add_ui(aux1, aux1, 1); 
	mpf_mul(aux1, aux1, y1);
	
    /* Calcula aux2 = aux1*2^n */
    mpf_set_ui(aux2, 2); 
	mpf_pow_ui(aux2, aux2, (2*k)+3); 
	mpf_mul(aux2, aux2, aux1);
	
    /* Calcula a1 = a0*(1+y1)^4 - aux2 */
    mpf_add_ui(a1, y1, 1); 
	mpf_pow_ui(a1, a1, 4); 
	mpf_mul(a1, a1, a0); 
	mpf_sub(a1, a1, aux2);

    mpf_clears(aux1, aux2, NULL);
}

/* Cálculo de Pi */
void calc_pi(mpf_t a, mpf_t pi)
{
    /* pi = 1/a */
    mpf_ui_div(pi, 1, a);
}

/* Borwein sequencial*/
void borwein(mpf_t pi, int *iter)
{
	mpf_t pi_ant;
    mpf_t a[2], y[2];
    int i;

    /* Inicialização do método */
	mpf_inits(a[1], y[1], pi_ant, NULL);
    start_values_borw(a[0], y[0]);

    /* Iteração para os cálculos */
    for(i=0; ;i++)
    {
        /* Calcula 'y1' e 'a1' */
        calc_y(y[1], y[0]);        
        calc_a(a[1], a[0], y[1], i);

        /* Armazena o pi anterior para convergência */
		mpf_set(pi_ant, pi);
		
		/* Calcula pi */
        calc_pi(a[1], pi);
	
		gmp_printf("Pi: %.*Ff\n", 15, pi);

        /* Verifica se convergiu */
        if(converge(pi_ant, pi)) 
			break;
			
		/* Armazena 'a' e 'y' para a próxima iteração */
		mpf_set(y[0], y[1]);
		mpf_set(a[0], a[1]);
    }

    mpf_clears(a[0], a[1], y[0], y[1], pi_ant, NULL);

    /* Número de iterações */
    *iter = i; 

}

/* Inicializa a struct para os cálculos paralelos */
void start_struct_borw(borw *aux)
{
	int i;

	aux->conv = FALSE;

    /* Primeiros valores do método */
    start_values_borw(aux->a[0], aux->y[0]);

	/* Vetor para armazenar todos os valores */
	for(i=1; i<16; i++)
    	mpf_inits(aux->a[i], aux->y[i], NULL);

	/* Inicializa as duas váriaveis que armazenam pi */	
	mpf_inits(aux->pi[0], aux->pi[1], NULL);
}

/* Limpa a struct */
void clear_struct_borw(borw *aux)
{
	int i;

    for(i=0; i<16; i++)
    	mpf_clears(aux->a[i], aux->y[i], NULL);

	mpf_clears(aux->pi[0], aux->pi[1], NULL);
}

/* Thread para cálculo de y */
void *thread_calc_y(void *a)
{
    borw *aux = (borw *) a;
	
    int id0 = (aux->i);	
	int id1 = (aux->i)+1;

    calc_y(aux->y[id1], aux->y[id0]);

    return NULL;
}

/* Thread para cálculo de a */
void *thread_calc_a(void *a)
{
    borw *aux = (borw *) a;

    /* Começa a armazenar 'a' a partir da segunda iteração */
    if(aux->i < 1) 
		return NULL;

    int id0 = (aux->i)-1;	
	int id1 = (aux->i);

    calc_a(aux->a[id1], aux->a[id0], aux->y[id1], aux->i-1);

    return NULL;
}

/* Thread para cálculo de Pi */
void *thread_calc_pi(void *a)
{
    borw *aux = (borw *) a;

    /* Começa o cálculo de pi a partir da terceira iteração */
    if(aux->i < 2) 
		return NULL;

	int id1 = (aux->i)-2;
	
	/* Armazena o pi anterior para convergência */
	mpf_set(aux->pi[0], aux->pi[1]);
	
	/* Cálculo de pi */
    calc_pi(aux->a[id1], aux->pi[1]);

    /* Se convergiu, a flag de convergência se torna verdadeira */
    if(converge(aux->pi[0], aux->pi[1])) 
		aux->conv = TRUE;

    return NULL;
}


/* Borwein paralelo */
void borwein_paral(mpf_t pi, int *iter)
{
    int i, j;
    borw b;
	pthread_t thread[3];

    start_struct_borw(&b);

    /* Realiza as iterações até atingir a precisão */
	for(i = 0; !b.conv; i++)
    {
        b.i = i;

        /* Uma thread é criada para 'a', 'y' e pi */
        pthread_create(&thread[0], NULL, thread_calc_y, &b);
        pthread_create(&thread[1], NULL, thread_calc_a, &b);
        pthread_create(&thread[2], NULL, thread_calc_pi, &b);

        /* Sincroniza as threads após elas terminarem */
        for(j = 0; j < 3; j++)
            pthread_join(thread[j], NULL);

		gmp_printf("Pi: %.*Ff\n", 15, b.pi[1]);
    }

	/* Número de iterações */
    *iter = i; 
	
    clear_struct_borw(&b);
}
